
// standard C++ headers
#include <iostream>
#include <fstream>
#include <string>

#include "problem.h"

// PROBLEMTYPE values:
//	1 : 1 group: 1 and 2 dimensions
//	2 : 2 group: 1, 2 and 3 dimensions
//	3 : 2 group: With upscattering
//	4 : 2 Group: IAEA Benchmark problem
//  5 : 2 Group: 2-D BWR problem by Demaziere
//	6 : 1 group: 1 and 2 dimensions with Neumann-Robin BC
//  7 : 2 Group: 2-D Benchmark from ANL
//  8 : 2 Group: 3-D test problem
#define PROBLEMTYPE 1

#if (PROBLEMTYPE == 1)

// 1 group eigen problem: 1 and 2 dimensions
// symmetric generalized eigen problem ; self-adjoint solution
#include "problem/problem1.cc"

#elif (PROBLEMTYPE == 2)

// 2 energy groups; no upscattering but with downscattering; 
// the loss matrix is block lower diagonal; unsymmetric matrix
// might not even be SPD 
#include "problem/problem2.cc"

#elif (PROBLEMTYPE == 3)

// 2 energy groups; with both upscattering and downscattering; 
// the loss matrix is full; unsymmetric matrix
// is not SPD ;
#include "problem/problem3.cc"

#elif (PROBLEMTYPE == 4)

// the standard IAEA 2D benchmark problem
// use an external mesh and solve the 2 group diffusion problem
// with downscattering accurately
#include "problem/iaea2d_problem.cc"

#elif (PROBLEMTYPE == 5)

// the standard IAEA 2D benchmark problem
// use an external mesh and solve the 2 group diffusion problem
// with downscattering accurately
#include "problem/bwr_problem.cc"

#elif (PROBLEMTYPE == 6)

// 1 energy group homogenous problem with Neumann and Robin BCs
#include "problem/problem6.cc"

#elif (PROBLEMTYPE == 7)

// 2 energy group problem benchmark from ANL
#include "problem/2d_benchmark_anl.cc"

#elif (PROBLEMTYPE == 8)

// 2 energy group problem benchmark from ANL
#include "problem/problem8.cc"

#else

#error "Provide a test problem."

#endif



#undef __FUNCT__
#define __FUNCT__ "Diffusion"
Real Problem::Diffusion(unsigned int ig, unsigned int imat) const
{
	libmesh_assert(ig < n_egroups) ;
	libmesh_assert(imat < n_materials) ;
	
	return xsdiff[imat][ig] ;
}

#undef __FUNCT__
#define __FUNCT__ "Removal"
Real Problem::Removal(unsigned int ig, unsigned int imat) const
{
	libmesh_assert(ig < n_egroups) ;
	libmesh_assert(imat < n_materials) ;
	
	return xsrem[imat][ig] ;
}

#undef __FUNCT__
#define __FUNCT__ "NuFission"
Real Problem::NuFission(unsigned int ig, unsigned int imat) const
{	
	libmesh_assert(ig < n_egroups) ;
	libmesh_assert(imat < n_materials) ;
	
	return xsnufiss[imat][ig] ;
}

#undef __FUNCT__
#define __FUNCT__ "NuFission"
Real Problem::Fission(unsigned int ig, unsigned int imat) const
{	
	libmesh_assert(ig < n_egroups) ;
	libmesh_assert(imat < n_materials) ;
	
	return xsfiss[imat][ig] ;
}

#undef __FUNCT__
#define __FUNCT__ "Scattering"
Real Problem::Scattering(unsigned src_group, unsigned int dest_group, unsigned int imat) const
{
	libmesh_assert(dest_group < n_egroups) ;
	libmesh_assert(src_group < n_egroups) ;
	libmesh_assert(imat < n_materials) ;
	
	return xsscatt[imat][dest_group][src_group] ;
}

#undef __FUNCT__
#define __FUNCT__ "EnergyPerFission"
Real Problem::EnergyPerFission (unsigned int ig) const
{	
	libmesh_assert(ig < n_egroups) ;
	
	return 1 ;
}

#undef __FUNCT__
#define __FUNCT__ "PromptFissionSpectrum"
Real Problem::PromptFissionSpectrum (unsigned int ig) const
{
	libmesh_assert(ig < n_egroups) ;
	
	return chip[ig]  ;
}

